library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE
library(mapview)
library(here)
## here() starts at C:/Users/jonah/Documents/GitHub/cw_rgis
# Read and export vector data ---------------------------------------------


# read a shapefile (e.g., ESRI Shapefile format)
# `quiet = TRUE` just for cleaner output
# for .shp, no need to call .dbf etc.
(sf_nc_county <- st_read(dsn = here("data/nc.shp"),
                         quiet = TRUE))
## Simple feature collection with 100 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
## Geodetic CRS:  WGS 84
## First 10 features:
##         county                       geometry
## 1         ashe MULTIPOLYGON (((-81.47258 3...
## 2    alleghany MULTIPOLYGON (((-81.23971 3...
## 3        surry MULTIPOLYGON (((-80.45614 3...
## 4    currituck MULTIPOLYGON (((-76.00863 3...
## 5  northampton MULTIPOLYGON (((-77.21736 3...
## 6     hertford MULTIPOLYGON (((-76.74474 3...
## 7       camden MULTIPOLYGON (((-76.00863 3...
## 8        gates MULTIPOLYGON (((-76.56218 3...
## 9       warren MULTIPOLYGON (((-78.30849 3...
## 10      stokes MULTIPOLYGON (((-80.02545 3...
# write to .shp and .gpkg
# 'append = FALSE' means the file will overwrite any previous file with that
# name, rather than add data points to it.
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.shp"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `ESRI Shapefile'
## Writing layer `sf_nc_county' to data source 
##   `C:/Users/jonah/Documents/GitHub/cw_rgis/data/sf_nc_county.shp' using driver `ESRI Shapefile'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
st_write(sf_nc_county,
         dsn = here("data/sf_nc_county.gpkg"),
         append = FALSE)
## Deleting layer `sf_nc_county' using driver `GPKG'
## Writing layer `sf_nc_county' to data source 
##   `C:/Users/jonah/Documents/GitHub/cw_rgis/data/sf_nc_county.gpkg' using driver `GPKG'
## Writing 100 features with 1 fields and geometry type Multi Polygon.
# save as .rds (not compatible with non-R software)
saveRDS(sf_nc_county,
        file = here("data/sf_nc_county.rds"))

# read a saved .rds file
sf_nc_county <- readRDS(file = here("data/sf_nc_county.rds"))


# Points ------------------------------------------------------------------

# This file contains a set of point data
(sf_site <- readRDS(here::here("data/sf_finsync_nc.rds")))
## Simple feature collection with 122 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.02452 ymin: 34.53497 xmax: -76.74611 ymax: 36.54083
## Geodetic CRS:  WGS 84
## # A tibble: 122 × 2
##    site_id                          geometry
##  * <chr>                         <POINT [°]>
##  1 finsync_nrs_nc-10013 (-81.51025 36.11188)
##  2 finsync_nrs_nc-10014 (-80.35989 35.87616)
##  3 finsync_nrs_nc-10020 (-81.74472 35.64379)
##  4 finsync_nrs_nc-10023  (-82.77898 35.6822)
##  5 finsync_nrs_nc-10024 (-77.75384 35.38553)
##  6 finsync_nrs_nc-10027 (-83.69678 35.02467)
##  7 finsync_nrs_nc-10029 (-80.65668 35.16119)
##  8 finsync_nrs_nc-10034 (-82.04497 36.09917)
##  9 finsync_nrs_nc-10041 (-80.46558 35.04365)
## 10 finsync_nrs_nc-10049 (-77.96322 34.58249)
## # ℹ 112 more rows
# View points
mapview(sf_site,
        col.regions = "black", 
        legend = FALSE)
# subsetting works the same as in any other data frame
(sf_site_f10 <- sf_site %>% 
    slice(1:10))
## Simple feature collection with 10 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -83.69678 ymin: 34.58249 xmax: -77.75384 ymax: 36.11188
## Geodetic CRS:  WGS 84
## # A tibble: 10 × 2
##    site_id                          geometry
##    <chr>                         <POINT [°]>
##  1 finsync_nrs_nc-10013 (-81.51025 36.11188)
##  2 finsync_nrs_nc-10014 (-80.35989 35.87616)
##  3 finsync_nrs_nc-10020 (-81.74472 35.64379)
##  4 finsync_nrs_nc-10023  (-82.77898 35.6822)
##  5 finsync_nrs_nc-10024 (-77.75384 35.38553)
##  6 finsync_nrs_nc-10027 (-83.69678 35.02467)
##  7 finsync_nrs_nc-10029 (-80.65668 35.16119)
##  8 finsync_nrs_nc-10034 (-82.04497 36.09917)
##  9 finsync_nrs_nc-10041 (-80.46558 35.04365)
## 10 finsync_nrs_nc-10049 (-77.96322 34.58249)
mapview(sf_site_f10,
        col.regions = "black",
        legend = FALSE)
# Lines -------------------------------------------------------------------

# This file contains line data (streams of Guilford County)
(sf_str <- readRDS(here::here("data/sf_stream_gi.rds")))
## Simple feature collection with 5012 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -80.04671 ymin: 35.9001 xmax: -79.53284 ymax: 36.25706
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-79.90748 36.18... fid000001
## 2  LINESTRING (-79.89768 36.20... fid000002
## 3  LINESTRING (-79.94756 36.15... fid000003
## 4  LINESTRING (-79.94152 36.25... fid000004
## 5  LINESTRING (-79.94206 36.20... fid000005
## 6  LINESTRING (-79.90003 36.18... fid000006
## 7  LINESTRING (-79.94737 36.20... fid000007
## 8  LINESTRING (-79.87595 36.18... fid000008
## 9  LINESTRING (-79.88022 36.25... fid000009
## 10 LINESTRING (-79.94859 36.25... fid000010
mapview(sf_str,
        color = "steelblue",
        legend = FALSE)
# The first ten segments are a random bunch of little lines
(sf_str_f10 <- sf_str %>% 
    slice(1:10))
## Simple feature collection with 10 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -79.94859 ymin: 36.15627 xmax: -79.87098 ymax: 36.25379
## Geodetic CRS:  WGS 84
##          fid                       geometry
## 1  fid000001 LINESTRING (-79.90748 36.18...
## 2  fid000002 LINESTRING (-79.89768 36.20...
## 3  fid000003 LINESTRING (-79.94756 36.15...
## 4  fid000004 LINESTRING (-79.94152 36.25...
## 5  fid000005 LINESTRING (-79.94206 36.20...
## 6  fid000006 LINESTRING (-79.90003 36.18...
## 7  fid000007 LINESTRING (-79.94737 36.20...
## 8  fid000008 LINESTRING (-79.87595 36.18...
## 9  fid000009 LINESTRING (-79.88022 36.25...
## 10 fid000010 LINESTRING (-79.94859 36.25...
mapview(sf_str_f10,
        color = "steelblue",
        legend = FALSE)
# Polygons ----------------------------------------------------------------

# Viewing the county data from earlier
mapview(sf_nc_county,
        color = "steelblue",
        col.regions = "tomato",
        legend = FALSE)
# Just Guilford County, using filter()
(sf_nc_gi <- sf_nc_county %>% 
    filter(county == "guilford"))
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -80.04238 ymin: 35.89106 xmax: -79.53034 ymax: 36.25032
## Geodetic CRS:  WGS 84
##     county                       geometry
## 1 guilford MULTIPOLYGON (((-79.53758 3...
mapview(sf_nc_gi)
# Mapping vector data -----------------------------------------------------

# You can put a map on a plot! And overlay multiple data sets
ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_site,
          color = "gray") +
  geom_sf(data = sf_str)

# but two layers covering the same area may not match up exactly...
ggplot() +
  geom_sf(data = sf_nc_gi,
          fill = "lightgreen")+
  geom_sf(data = sf_str,
          color = "steelblue")

# Exercise ----------------------------------------------------------------

# Read stream line data for Ashe County
sf_str_as <- readRDS(file = here("data/sf_stream_as.rds"))

# Confirm that the CRS matches the county data set
print(sf_str_as) # CRS is WGS 84
## Simple feature collection with 6361 features and 1 field
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -81.73891 ymin: 36.24445 xmax: -81.24562 ymax: 36.588
## Geodetic CRS:  WGS 84
## First 10 features:
##                          geometry       fid
## 1  LINESTRING (-81.32158 36.47... fid000001
## 2  LINESTRING (-81.31133 36.47... fid000002
## 3  LINESTRING (-81.38594 36.45... fid000003
## 4  LINESTRING (-81.39307 36.46... fid000004
## 5  LINESTRING (-81.3454 36.533... fid000005
## 6  LINESTRING (-81.34719 36.53... fid000006
## 7  LINESTRING (-81.3328 36.518... fid000007
## 8  LINESTRING (-81.34964 36.50... fid000008
## 9  LINESTRING (-81.35949 36.53... fid000009
## 10 LINESTRING (-81.35094 36.52... fid000010
print(sf_nc_county) # CRS is WGS 84
## Simple feature collection with 100 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.45662 ymax: 36.58973
## Geodetic CRS:  WGS 84
## First 10 features:
##         county                       geometry
## 1         ashe MULTIPOLYGON (((-81.47258 3...
## 2    alleghany MULTIPOLYGON (((-81.23971 3...
## 3        surry MULTIPOLYGON (((-80.45614 3...
## 4    currituck MULTIPOLYGON (((-76.00863 3...
## 5  northampton MULTIPOLYGON (((-77.21736 3...
## 6     hertford MULTIPOLYGON (((-76.74474 3...
## 7       camden MULTIPOLYGON (((-76.00863 3...
## 8        gates MULTIPOLYGON (((-76.56218 3...
## 9       warren MULTIPOLYGON (((-78.30849 3...
## 10      stokes MULTIPOLYGON (((-80.02545 3...
# Export as .gpkg
#st_write(sf_str_as,
#         dsn = "data/sf_stream_ashe.gpkg",
#         append = FALSE)
          # error: wrong field type for fid

# Map streams along with county boundaries
ggplot() +
  geom_sf(data = sf_nc_county) +
  geom_sf(data = sf_str_as)

# Subset county layer to Ashe County and remap
sf_nc_as <- sf_nc_county %>% 
  filter(county == "ashe") # subsetting

ggplot() +
  geom_sf(data = sf_nc_as) +
  geom_sf(data = sf_str_as,
          color = "steelblue") # re-mapping